program coaxial
  use constants
  use select_precision
  use mymodule
  use dispmodule
  implicit none

  real(wp)                                :: a,b,c,d,alpha,p,s,delta,epsilon_r,sigma,freq,delta_s,omega,Rs,dZ0
  real(wp),dimension(:,:),allocatable     :: stats,  res_vector
  real(wp),dimension(:),allocatable     ::  d_inner , d_outer

  integer(ip)                             :: N,i, j,Ngeo
  real(wp),dimension(4)                     :: answer
  integer(ip),parameter                     :: file_handle = 100
  
  open( unit = file_handle , file='coaxial.out' , action='write' )
  
  !!!!!!!!!!!!!!!
  ! CONSTANTS
  !!!!!!!!!!!!!!!
  a = .1_wp
  b = .1_wp
  c = 1._wp
  d = 1._wp
  p = 5.e-3_wp
  s = 5.e-3_wp
  alpha =  30._wp*DEG2RAD
  epsilon_r = 1._wp
  delta = c/10._wp
  freq = 1.0e9_wp
  sigma = 5.813e7_wp
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Skin depth calculation
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  omega = 2._wp*PI*freq
  delta_s = sqrt( 2._wp/( omega*MU0*sigma ) )
  Rs = 1._wp/( delta_s*sigma )
    
  print *,'Give number of  times to evaluate impedance:'
  read *,N
  if (N<3) then 
    N = 3
  endif
  
  !print *,'Give number of rows in table:'
  !read *,Ngeo
  Ngeo = 2
  
  allocate( d_outer(Ngeo)  , d_inner(Ngeo) , res_vector(Ngeo,3) )
  allocate( stats(N,4)  )
  
  res_vector = ZERO
  call linspace( d_inner , 2.e-3_wp , 2.e-3_wp - delta_s*5._wp )
  call linspace( d_outer , 4.0e-3_wp , 4.0e-3_wp + delta_s *5._wp )
  delta = d_outer(1)/10._wp
  !call disp( res_vector)
  stats = ZERO
  
  do j = 1,Ngeo
    stats = ZERO

    do i=1,N
        stats(i,2) = delta/TWO**(i-1)
        answer    = impedance( d_inner(j),d_inner(j) , d_outer(j), d_outer(j) , p,s,alpha,epsilon_r, delta/TWO**(i-1) )
        stats(i,1) = answer(1)
        !print "(a,i5,a,f9.5,a,2f12.3)"  , 'Solving ' ,int(answer(2),ip) ,' unknowns. Impedance is' , answer(1) &
        !       , ' Ohms. Fill/solve:' , answer(3),answer(4)
    enddo
  
    do i=3,N
        stats( i , 3 ) = log( (stats(i-2,1)-stats(i-1,1))/(stats(i-1,1)-stats(i,1)) )/log(TWO)
        stats( i , 4 ) = stats(i,1) - ( stats(i-1,1) - stats(i,1) )/(TWO**stats(i,3) - ONE )
    enddo
    !call disp( stats ,fmt='f16.9')
    res_vector( j , 1 ) = d_inner(j)
    res_vector( j , 2 ) = d_outer(j)
    res_vector( j , 3 ) = stats( N , 4 )
  enddo
  
  dZ0 =  ( res_vector(2,3) - res_vector(1,3) )/10._wp
  print "(a,f6.4,a)", 'Loss alpha_c: ' , Rs/ETA0/delta_s*dZ0/res_vector(1,3)*NPR2DB,' [dB/m]'
    
  a = d_inner(1)/2._wp
  b = d_outer(1)/2._wp
  print  "(a,f6.4,a)", 'Analytic alpha_c: ' , Rs/(2._wp*ETA0*log(b/a))*(1/b+1/a)*NPR2DB, ' [dB/m]' 
  
  
  !call disp( res_vector ,fmt = 'f25.16' )

  call disp( res_vector , unit=file_handle , fmt = 'f25.16' )
  call disp( res_vector , fmt = 'f25.16' )
  
  close( file_handle )
  deallocate( stats , res_vector , d_outer , d_inner)  
  end program coaxial
  
  
